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Abstract. We study the effects of periodically time varying mass on the stability of the Helmholtz oscillator, 
which, when linearised, takes the form of Ince's equation and exhibits parametric resonance. The resonance 
f"""""" ' regions in the parameter space are mapped and we use the Melnikov function to demonstrate that the 

C D ' parametric instability does not affect the qualitative dynamics of the nonlinear system. 

o : 

1. INTRODUCTION 

C . 

C3 ■ The recent activity in the application of Nonlinear Dynamics to problems in Naval Architecture is exemplified 
' 4 by the papers in [13j . One particular success has been a more rational understanding of beam sea capsize 
. and the development of a robust mathematical toolkit for the study of this problem. 

Following the work of Wright and Marshfield, [20J it has become common to study the roll dynamics using 
. a one degree of freedom nonlinear oscillator taking the form 
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(1) 6 + g(6)+G z (8) = F(t). 



Here 9 is the roll angle relative to the wave normal, g(9) is the, often linear- viscous, damping function, and 
Gz{0) (so denoted as it is usually some polynomial fit to the GZ curve) encapsulates the restoring moment. 
' The forcing, F(t), is often taken to be a simple sinusoidal function, see for example [17], although stochastic 
Q . models of forcing are also studied, e.g. [3]. 

' Observing that capsize is equivalent to the roll angle exceeding the angle of vanishing stability leads to the 
_^ . consideration of the more general dynamical problem of escape of a particle from a potential well. The 
J> ' Helmholtz oscillator 

(2) x + fix + x - x 2 = Fsincot, 



\ known to Naval- Architects as the Helmholtz -Thompson equation, provides a simple archetype to study 



the phenomena of escape directly. It has been central to a number of developments. Thompson's paper 
! |14j provides an overview of the dynamic response of the system which has been studied experimentally by 
c"^i Gottwald et. al. [2 . The inhibition of chaotic escape is considered in the context of © in [T] and a more 
general approach to the problem is discussed in [BJ. The equation finds direct application in the study of 
bubble dynamics [5] and is much discussed in the Naval Architecture literature, see [TBJ. The engineering 
integrity diagram, and the use of Melnikov theory to predict parameter values for which erosion of basins 
of attraction takes place were developed in this context. These concepts continue to find fruitful application 
in quantification of capsize resistance, see [T2"] . 



Providing a reasonable nonlinear dynamical theory of the stability of a damaged ship remains a considerable 
^ • challenge. Compared with intact stability there are a number of dynamical differences which effect ship 
motion. The mass of a progressively flooding ship changes with time. The motion becomes coupled with 
that of water on deck or in compartments. The water itself provides non-trivial modelling problems. Masses 
of water lead to biases which may also vary with time. 

The increased complexity means that one is forced to make a choice between sacrificing dynamical information 
by making simplifying assumptions or dynamical insight by building less simple comprehensive models. 
When making assumptions to simplify a model enough that some theoretical analysis is possible there is a 
danger of assuming so much that all physical meaning is lost. The balance between finding a model that is 
simple enough to analyze yet still manages to capture the essence of some real physical dynamics is difficult 
to attain. One successful example is found in Murashigc et al's analysis of flooded ship motion, see [5J and 
references therein. 
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In this paper we choose a different approach. The dynamic phenomena listed above will have varying 
importance in building a comprehensive theory of damage stability. Rather than start with a physical model 
and make simplifying assumptions, we have chosen to put one of these effects, time varying mass, directly 
under the microscope. Whilst some properties of time varying mass dynamics which will influence the 
behaviour of any model whose mass has a periodic component are discussed, we make no claim to direct 
physical relevance. 

The model we study is a Helmholtz oscillator with time varying mass. We use the generalised form of 
Newton's Second Law, 

F = dt 

(3) = m(t)x + m(t)x , 

so that an oscillator modelled as mi + g(x,x) = F(t) whose mass varies with time becomes 

(4) m(t)x + g(x,x) + m(t)x = F(t). 
Equation @ becomes 

, „ 8 + m(t) x — x 2 sinw/rt 

(5) x + - r + 1 x + — r^ = F rr~ 

v ' m(t) m(t) m(t) 

where we have divided through by m(t). 

In this paper we study mass varying sinusoidally with time and choose the simple form 

(6) m(t) = 5 + 7 sincj m £. 

Substituting © into ([5]) and writing 5 = 1/8 (in order to simplify the perturbation calculations which follow) 
we have 

.. , s/3 + 7W m cosw m t . j x-x 2 j sinu F t 

(7) x + o — x + o — = or — . 

1 + "/5 sin uj m t l+7^sino; m i 1 + 7$ smu) m t 

2. LINEAR STABILITY ANALYSIS OF FIXED POINTS 

Our first step is to consider the linear stability of the two fixed points of the system with F = 0. We start 
by scaling time in order to normalise the frequency of the mass variation. Setting r = ui m t yields 

(s\ ■■ , ; £+7cost . 5 x-x 2 

(8) x + d— x-\ : = 0. 

I + 76 sin t UJ m 1 + 7(5 sin t 

after relabelling t = t. Writing ([8|) as a system 
(9a) x = y 

5 x-x 2 - — + 7 cos t 

(9b) y = -^-^ 5^ 1 y 

1 + j 6 sin t 1 + 7<5 sin t 

we see that the two fixed points are (x, y) — (0, 0) and (1, 0). We know that for 7 = the origin is a stable 
focus or centre as (3 > or /? = respectively. The point (1,0) is a saddle independent of the value of (3. 
For both fixed points the linearised system takes the form 

x = A(t)x 

where A(t) = A(t + 2ir) is a periodic matrix and we are in the setting of Floquet theory [4]. For the saddle 
point we have 



(10) 4(1,0) (*) = 

The origin gives 




1+7(5 sin t l+7<5 sin t 







(11) 4(0,0) (*) = I I ^+7 cost 

\ 1+7^ sin t 1+7^ sin t 
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THE FIXED POINT AT (1,0). The linear stability analysis of this fixed point does not require the full 
machinery of Floquet theory. The ODE linearised about this point is 

x H x -. —x = v. 

l + 7(5sini (1 + 7 <5 sin £ 



Upon expanding x = x n + -fxi + 0(7 2 ) and 5 = <5 + ySi + 0{^ 2 ) and collecting terms in 7 we see that the 
O(l) equation is 

x H x a - OqXq = 0. 

The origin in this equation is stable only when Sq < which is physically meaningless. We conclude that the 
fixed point of the nonlinear system is unstable for all S > 0, regardless of the value of j3 and u> m . 

THE FIXED POINT AT THE ORIGIN. The ODE linearised about the origin is 

-f-+(57Cosi 5 

(12) x + ^ , x + ? ^x = 0. 

l + 7 <5sint w r 2 „ (1 +7^sint) 



and in this case there are a number of regions where a Mathieu type instability occurs. It follows from 
Floquet theory that transition curves in the parameter space correspond to solutions of the system which 
have periods equal to and double the fundamental period of the system. Thinking about periodicity is easier 
when the fundamental periods under consideration are 2ir and At: so we scale time again setting t — 2r. As 
we wish to follow Rand's treatment of Ince's equation 

(1 + a cos 2t)x + b sin 2tx + (c + d cos 2t)x = 0, 

see [TU], we also translate time r h- > r + 2ir. Relabelling r — t for notational convenience we have 

(13) (1 + 57 cos 2t)x + 2 (— - $7 sin 2t] x + A-^-x = 0. 



When — this is Ince's equation with a — 6~f, b = —25"/, c = A-Jr and d = 0. We will consider /3 ^ when 
the P = dynamics are understood. 

Remark. The scaling of time has no effect on the stability of the origin in the linearised system. The stable 
parameters for system ([T2"j) are exactly those that are stable for system (TTBl . 

As transition curves are characterised by the existence of solutions of period 2ir and Air we use harmonic 
balance to seek solutions of this form. Substituting x = X^^Lo a ™ cos n ^ b n sinnt into (|13|) . combining 
trigonometric terms and equating coefficients of cos nt and sin nt leads to an infinite set of linear homogeneous 
equations in the Fourier coefficients and bi. The first few of these are 



(57 



4— a = 



ai - 2<^7 a 3 = 



(57 6 \ 3~ 
( aAt- - A I a 2 + -A5ja 4 = 



The existence of a transition curve in the parameter space is equivalent to the existence of a nontrivial solution 
to this system of equations. This in turn is equivalent to the vanishing of the associated determinants of the 
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system. This last condition says that the transition curves in the parameter space are composed of points 
which satisfy the equations 



(14a) 



(14b) 



(14c) 
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We note that 
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A„ 2 , =4— A b2i . 

This fact has two consequences. First of all it is clear that 5 — is a solution of (|14a|) so that the first 
transition curve we find is 



(15) 



6 = 0. 



Secondly all other zeros of A Q2i = will coincide with the zeros of Af, 2 
that these curves emanate from the points 



0. Upon setting 7 = it is clear 



6 = 



2nu) r , 



e{i,2,3,...}. 



However because they coincide the regions they bound have zero width and no instability occurs. The 
coincidence of two transition curves and associated vanishing of a regime of instability is known as coexistence 
see [7] or [10] for details. 

Remark. The lack of instability is not a consequence of the vanishing width of the unstable regime. Periodic 
solutions exist on the transition curves and those on the A a2i = are even functions of t. On the Afc 2i = 
the solutions are odd. As odd and even solutions are linearly independent we have, on two coinciding curves, 
a spanning set for the vector space of solutions. A solution space spanned by bounded solutions cannot 
contain an unbounded solution. 
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Figure 1 . The first three regions of parametric instability. 



In order to evaluate the solution curves of A 



02i + l 



and A& 2i+1 



we substitute 



5 = S + 5x1 + $2j 2 + hi 3 + <^47 4 + ■ • • 

into a finite, n x n where n is a small integer, truncation of the determinants A Q2i+1 and Aj, 2i+1 , equate to 
zero and collect coefficients of 7 in the usual fashion. Using Maple and determinants of order 25 x 25 we can 
easily write down the first few curves. 

§7 + 0< 7 *) 



(16a) 


.5 = 


4 


(16b) 


5 = 


4 


(16c) 


5 = 


9 2 


(16d) 


5 = 


9 2 


(16e) 


5 = 


25 2 


(16f) 


5 = 


25 2 



16767 fi n 

4096 ml 
16767 fi n 

4096 ml 
1109375 



6561 



<7 3 



131072 

6561 o t, 

oj 8 ^ 3 

131072 mr 



0(7 4 ) 
0(7 4 ) 



^ 6 7 2 + 



3030048828125 



12288 ' 1811939328 
1109375 6 2 3030048828125 in 4 

■ W 7 + 1Q nnonooo W m 7 



12288 ' 1811939328 
As the equations of the curves suggest, it is possible to prove: 

Proposition 1. The pair of transition curves emanating from the point 

'(2k-l)w. 



2197265625 

2147483648 
2197265625 

2147483648 



^V + 0( 7 b ) 
^™ 2 7 5 + 0(7 6 ) 



coincide to order 7 



2fc-2 



They split at order 7 



2fe-l 



The geometry of this fact is that as 5 grows the curves becomes thinner and thinner. The areas of the 
parameter space in which instability occurs decreases. Figure[T]illustrates the first three regions of instability 
all of which appear thin on this scale. A closer look at the first region (along with the effect of (3 > 0, discussed 
below) is seen in Figure [U 

Remark. We have not mentioned that a mathematically rigourous treatment would have to address the 
divergence of the determinants in their current form. The solution to this problem involves noting that 
the individual coefficients of 7™ are defined as ratios of infinite determinants. One must then consider the 
coefficient as the ratio of two k x k determinants, let k — > 00, and show that the expression converges. 
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Figure 2. The First Region of Instability (With and Without Damping). u) m = 1 



This, along with the proof of Proposition [TJ are essentially exercises in algebra and we do not reproduce the 
calculations here. 

When (3 ^ the above analysis fails as the systems of equations for the coefficients of the Fourier series 
do not decouple in a nice way. We use a perturbation argument instead to illustrate the effect of /3. It is 
necessary to assume that (3 is of the same order as 7 which we consider to be a small parameter. Thus we 
write (3 = 7/3 and we are studying 



(17) 



x + jS 



-2- + cos t 
1 + 7<5 sini 



1 + 7<5 sin t 



0. 



We expand 5 = S + Sij + 5 2 j 2 + CK7 3 ), x = x + xij + x 2 j 2 + 0(j 3 ), substitute and collect coefficients in 
the usual way. Considering <5 = w^„/4 we see that 



£ 1 = ± 



32 



16/5 2 . 



Now remembering that 7/? = (3 we have 



(18) 



4 32 



16(3 2 + 0(Y)- 



In Figure [U we see that the effect of (3 > is to lift the curves from the 7 = axis, increasing the amplitude 
of mass variation that is required to produce an instability. Figure G3 illustrates the numerical response of 
the system in (x, t) space. The left hand graph corresponds to parameters outwith the region of parametric 
resonance, the right has parameters inside the region. 

3. THE MELNIKOV FUNCTION FOR THE NONLINEAR SYSTEM 

Having understood the linear dynamics of the two fixed points we now wish to consider the full nonlinear 
problem. We will use Melnikov Theory (see [T5] or [T!|]) to understand which parameter values support 
homoclinic tangling (and therefor safe basin erosion) and compare the analytic results with numerically 
calculated integrity curves. 

To calculate the Melnikov Function we need to consider ([7|) as a nearly integrable system taking the form 
(19) x = f(x) + eg(x,i). 

We assume that F, 7, and (3 are O(e) and expand 

5 



1+7(5 sin ui m t 



1 — ^8smuj m t + 0(7 ). 



STABILITY ANALYSIS OF THE HELMHOLTZ OSCILLATOR WITH TIME VARYING MASS. 



7 




5 hat = . 25, y =.038 



5 hat = .25, y=.042 x 



Figure 3. Trajectories of the Linear System. The left hand solution has parameters just 
outside the resonant region. The right hand, just inside. In both cases the solution were 
plotted using the Matlab ode45 routine with initial condition (1,0). 



Writing F — eF, 7 = ej and (3 = eft for notational convenience we have ([7]) in the form (fTT)]) and the Melnikov 
function is 



M(t ) 



FS sin (cj F (t + 1 )) y h (t) - 6/3yl(t) - &yw m cos (u) m (t + 1 )) y 2 h {t) 

+ j5 2 sin (Lu m {t + t ))(x h {t) - x 2 h (t))y h {t) \dt. 



Here 



Xh(t) 



1 



Vh = Xh 



1 + cosh 

is the parameterised orbit homoclinic to (1,0). Using trigonometric reduction, changes of variables and 
residue calculus this evaluates as 



(20) 



M(t ) 



QFljJpTT COS UJFto 

sinh^f 



-5.0. 



7 



3 7tcj4 coso; m to 
5 5 sinh ^ 



Notice that when 8=1 and 7 = this is identical to the Melnikov function for the standard Hclmholtz 
oscillator, [T5] . 

Remark. In cases where u> m ^ lu f and in particular when Lu m /u>F ^ Q we must rely on the generalisation of 
Melnikov theory to quasi-periodic perturbations due to Meyer and Sell [S] and Wiggins [18 . In this paper 
Melnikov Function is used as a synonym for Generalised Quasi-Periodic Melnikov Function. 

The Effect of Linear Resonance on Basin Erosion. To estimate the minimal F for which erosion is expected 
we set M(to) = and solve the function for F, selecting to to minimise as appropriate. The first question 
we wish to answer is does basin erosion differ depending on whither or not we are located within a region 
of parametric instability. Consider the region illustrated in Figure [H In this case uj m = 1 and the region is 
centred at 5 = 0.25. First we fix j3 = 0.01, S — 0.25 and vary 7 to cross the boundary of the region. With 
these parameters the Melnikov estimate of F looks like 

(.0015 sinh 2tt - 2.25^) sinh 2ttu f 



Fm{i,u f ) 



6ttu; f sinh 27r 



and for fixed iv F this is just a linear function of 7, see Figure HI 

If the presence of parametric instability had any interesting effect upon basin erosion we would expect some 
effect on Fm as 7 moved across the stability boundary. There is, however, no cusp, blip, change in gradient 
or other effect in the Melnikov estimate as the parameters move across the transition curve. 
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eo F = 1.15 
COp = 1 

CQp = 0.85 



-5 

0.03 



0.045 

Y 



0.05 0.055 



Figure 4. Fm as a function of 7. Even as 7 passes over the stability boundary at 7 = 0.04 
the function remains linear. 

The parametric resonance has no effect on the speed of homoclinic tangling and therefor no import for basin 
erosion. As the tangling and erosion processes are governed by the invariant manifolds of the saddle this 
result is unsurprising. Intuitively the homoclinic orbit does not see the effect of the local dynamics at the 
other fixed point. Time varying mass is, however, important in its quantitative effects on basin erosion. 

Time Varying Mass as Excitation. Notice in Figure |J] that when 7 is slightly greater than 0.055 the lines all 
cross Fm = 0. To understand this we consider the Melnikov Function with F = as an estimator of 7. We 
have 

3 TTLO^ COS LO m t 



M(* ) 



■-5*0- 
5 



7 



5 $ sinh ■ 



(< 2 



The energy balance between the linear damping and the excitation provided by mass variation which yields 
a zero of the Melnikov Function occurs when 



(21) 



2P/3sinh^ 
vs 



u^- 8 2 ) cos uj m t 



For the parameters in Figure[4]this evaluates to 7 ~ .05681, so for 7 above this value the mass variation alone 
is enough to cause basin erosion. It is worth noting that this also occurs away from regions of parametric 
instability. For uj m = tt there is no parametric resonance near 6=1, yet when F = 0, 7m ~ 6.467/3 still 
provides a zero of M(0). 



The General Effect of Time Varying Mass. The form of the Melnikov Function for periodic time varying 
mass (|20[) can be considered as the usual Helmholtz-Thompson Melnikov function plus a correction which 
takes into account the mass variation. The mass varying term is 



"7 



3 ttuj 1 ^ COS UJ m t 
5~f~ sinh ^sT 



and this changes sign at 6 = ±w„. This means that when S = uj^ the time varying mass does not affect the 
Melnikov estimate of the first homoclinic tangle. Figure [5] is the engineering integrity diagram for this case. 
Although there is some minor quantitative effect on the size of the initial safe basin due to 7 > the position 
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of the dover cliff reduction in basin area is not effected. If we set M(to) = and solve for F as usual we get 

sinh^ o sinh^ 

(22) F= — I 7 ^ ^ cos ^ m t 

57tcj^ cosw^to 107rw?, sinh coso^ t \ "v ' 

Assuming 7 > this makes it clear that when <5 2 > cj„ the effect of the time varying mass is to delay basin 
erosion. When S 2 < the basin erosion occurs for smaller F. This effect is obvious in the integrity diagrams 
in Figures [H and [7] 



4. CONCLUSIONS 

When Time Varying Mass in an oscillator is taken to be periodic the effect is similar to parametric forcing. 
Mathieu type resonance occurs although due to coexistence only half as many regimes of parametric resonance 
exist. Also as the parameter 5 grows, i.e. as the mean mass of the oscillator tends to zero, the regimes 
of instability become thinner and thinner. Whilst the linear parametric resonance does not have a direct 
qualitative effect on the dynamics of basin erosion, the quantitative effects are significant and depend critically 
on the relationship between the mean mass and the frequency of the variation. In some cases it can be broadly 
considered a stabilising influence whilst in others the opposite is true. 

In a realistic model of a damaged ship the assumption of (sinusoidal) periodicity for the variation of the mass 
due to water ingress is unlikely to be realised exactly. However it is often the case that flood water mass has 
some periodic component. When this occurs the effects studied in this paper are likely to be an important 
element of a complete picture of the dynamics. 
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Figure 5. Integrity Diagram for 6 = uif n = 1. The area of the safe basin (i.e. the union 
of the basins of attraction of the systems bounded attractors) normalised with respect to 
F = 7 = is plotted against the forcing parameter. 
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Figure 6. 8 = l,w m = 0.85 Figure 7. 5 = l,w m = 1.15 
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